# coding=utf-8
#*****************************************************************************
#       Copyright (C) 2008 Franco Saliola <saliola@gmail.com>,
#                          Sébastien Labbé <slabqc@gmail.com>,
#                          Arnaud Bergeron <abergeron@gmail.com>,
#                          Amy Glen <amy.glen@gmail.com>
#
#  Distributed under the terms of the GNU General Public License version 2 (GPLv2)
#
#  The full text of the GPLv2 is available at:
#
#                  http://www.gnu.org/licenses/
#*****************************************************************************
from itertools import cycle, count
from random import randint
from sage.structure.sage_object import SageObject
from sage.words.word import Words, FiniteWord, is_Word, is_FiniteWord
from sage.words.word_content import MakeContent
from sage.words.alpha import Alphabet
from sage.words.utils import *

w_0_1 = Words([0, 1])

def _build_tab(sym, tab, W):
    r"""
    Internal function building a coding table for the phi_inv_tab function.
    
    TESTS:
        sage: from sage.words.word_generators import _build_tab
        sage: _build_tab(1, [], Words([1, 2]))
        [1]
        sage: _build_tab(1, [1], Words([1, 2]))
        [1, 2]
        sage: _build_tab(2, [1], Words([1, 2]))
        [2, 2]
        sage: _build_tab(2, [1, 2], Words([1, 2]))
        [2, 2, 1]
        sage: _build_tab(1, [2, 2], Words([1, 2]))
        [1, 1, 2]
    """
    res = [sym]
    if len(tab) == 0:
        return res
    if sym == 1:
        res += tab
        res[1] = (res[1] % len(W.alphabet())) + 1
        return res
    w = W([sym]).delta_inv(W, tab[0])
    w = w[1:]
    res.append((w[-1] % len(W.alphabet())) + 1)
    for i in xrange(1, len(tab)):
        w = w.delta_inv(W, tab[i])
        res.append((w[-1] % len(W.alphabet())) + 1)
    return res

def _nth_digit(n):
    r"""
    Returns the digit in position n of the Thue-Morse word.
    
    The nth digit of the Thue-Morse word can be defined as the number
    of bits in the 2-complement representation of the position
    modulo 2 which is what this function uses.  The running time
    is O(lg n) where n is the position desired.
    
    INPUT:
        n -- integer, the position
        
    OUTPUT:
        0 or 1 -- the digit at the position
        
    TESTS:
        sage: from sage.words.word_generators import _nth_digit
        sage: _nth_digit(0)
        0
        sage: _nth_digit(3)
        0
        sage: _nth_digit(32)
        1
    """
    for tn in count():
        if n == 0:
            return tn & 1
        n &= n - 1


class WordGenerator(object):
    
    def ThueMorseWord(self):
        r"""
        Returns the Thue-Morse word t on the letters 0 and 1.

        There are several ways to define t. We use the following definition:
        t[n] is the sum modulo 2 of the digits in the binary expansion of n.

        EXAMPLES:
            sage: t = words.ThueMorseWord(); t
            Thue-Morse word on the alphabet [0, 1]
            sage: print t[:50]
            word: 01101001100101101001011001101001100101100110100101
        """
        w = w_0_1(_nth_digit, format="function")
        w.rename("Thue-Morse word on the alphabet %s" % w.alphabet().string_rep())
        return w


    def ChristoffelWord(self,p,q):
        r"""
        Returns the lower Christoffel word of slope p/q, where p and q are
        relatively prime nonnegative integers.

        The Christoffel word of slope p/q is obtained from the Cayley graph of
        Z/(p+q)Z with generator q as follows. If u --> v is an edge in the
        Cayley graph, then v = u+p mod (p+q). Label the edge u --> v by 1 if u
        < v and 0 otherwise. The Christoffel word is the word obtained by
        reading the edge labels along the cycle beginning from 0.

        EXAMPLES:
            sage: w = words.ChristoffelWord(1,0); w
            Lower Christoffel word of slope 1/0 over the alphabet [0, 1]
            sage: print w
            word: 1
    
            sage: w = words.ChristoffelWord(0,1); w
            Lower Christoffel word of slope 0/1 over the alphabet [0, 1]
            sage: print w
            word: 0
    
            sage: w = words.ChristoffelWord(1,1); w
            Lower Christoffel word of slope 1/1 over the alphabet [0, 1]
            sage: print w
            word: 01
    
            sage: w = words.ChristoffelWord(4,7); w
            Lower Christoffel word of slope 4/7 over the alphabet [0, 1]
            sage: print w
            word: 00100100101
        """
        return ChristoffelWord_Lower(p, q)
        
    LowerChristoffelWord = ChristoffelWord

    def UpperChristoffelWord(self,p,q):
        r"""
        Returns the upper Christoffel word of slope p/q, where p and q are
        relatively prime nonnegative integers.

        The upper Christoffel word of slope p/q is equal to the reversal of the
        lower Christoffel word of slope p/q. Equivalently, if xuy is the lower
        Christoffel word of slope p/q, then yux is the upper Christoffel word
        of slope p/q (because u is a palindrome).

        EXAMPLES:
            sage: w = words.UpperChristoffelWord(1,0); w
            Upper Christoffel word of slope 1/0 over the alphabet [0, 1]
            sage: print w
            word: 1
    
            sage: w = words.UpperChristoffelWord(0,1); w
            Upper Christoffel word of slope 0/1 over the alphabet [0, 1]
            sage: print w
            word: 0
    
            sage: w = words.UpperChristoffelWord(1,1); w
            Upper Christoffel word of slope 1/1 over the alphabet [0, 1]
            sage: print w
            word: 10
    
            sage: w = words.UpperChristoffelWord(4,7); w
            Upper Christoffel word of slope 4/7 over the alphabet [0, 1]
            sage: print w
            word: 10100100100
        """
        w = words.LowerChristoffelWord(p, q).reversal()
        w.rename("Upper Christoffel word of slope %s/%s over the alphabet %s" % (p, q, w.alphabet().string_rep()))
        return w

    def FibonacciWord(self, construction_method="recursive"):
        r"""
        Returns the Fibonacci word on the letters 0 and 1.

        INPUT:
            construction_method -- can be any of the following:
                    "recursive", "fixed point" (see below for definitions).

        Recursive construction: the Fibonacci word is the limit of the
        following sequence of words: S_0 = 0, S_1 = 01, S_n = S_{n-1} S_{n-2}
        for n >= 2.
        
        Fixed point construction: the Fibonacci word is the fixed point of the
        morphism 0 -> 01 and 1 -> 0. Hence, it can be constructed by the
        following read-write process: 

            beginning at the first letter of 01, 
            if the next letter is 0, append 01 to the word;
            if the next letter is 1, append 1 to the word;
            move to the next letter of the word.

        EXAMPLES:
            sage: from sage.words.word_generators import w_0_1
            sage: from math import floor, sqrt
            sage: w = words.FibonacciWord(construction_method="recursive"); w
            Fibonacci word over [0, 1], defined recursively
            sage: print w[:99]
            word: 010010100100101001010010010100100101001010010010100101001001010010010100101001001010010010100101001
            sage: u = words.FibonacciWord(construction_method="fixed_point"); u
            Fibonacci word over [0, 1], defined as the fixed point of a morphism
            sage: print u[:99]
            word: 010010100100101001010010010100100101001010010010100101001001010010010100101001001010010010100101001
            sage: golden_ratio = (1 + sqrt(5))/2.0
            sage: a = golden_ratio / (1  + 2*golden_ratio)
            sage: wn = lambda n : int(floor(a*(n+2)) - floor(a*(n+1)))
            sage: f = w_0_1(wn); f
            Infinite word over [0, 1]
            sage: f[:10000] == w[:10000]
            True
            sage: f[:10000] == u[:10000]
            True
        """
        if construction_method == "recursive":
            def RecursiveConstructionIterator():
                Fib0 = [0]
                Fib1 = [0,1]
                n = 0
                while True:
                    it = iter(Fib1[n:])
                    for i in it:
                        n += 1
                        yield i
                    else:
                        Fib1, Fib0 = Fib1 + Fib0, Fib1
            w = w_0_1(RecursiveConstructionIterator(), format='iterator')
            w.rename("Fibonacci word over %s, defined recursively" % w.alphabet().string_rep())
            return w
        elif construction_method == "fixed_point":
            w = self.FixedPointOfMorphism({1:[0],0:[0,1]}, 0)
            w.rename("Fibonacci word over %s, defined as the fixed point of a morphism" % w.alphabet().string_rep())
            return w
        else:
            raise NotImplementedError

    def FixedPointOfMorphism(self, morphism, first_letter):
        r"""
        Returns the infinite word beginning with first_letter that is the fixed
        point of the given morphism. Here morphism is a dictionary.

        Recall that the fixed point w of morphism is defined as
        
            morphism(w[0]) morphism(w[1]) morphism(w[2]) ...

        EXAMPLES:
            sage: mu = {0:[0,1], 1:[1,0]}
            sage: tm = words.FixedPointOfMorphism(mu,0); tm
            Fixed point beginning with 0 of the morphism {0: [0, 1], 1: [1, 0]}
            sage: TM = words.ThueMorseWord()
            sage: tm[:10000] == TM[:10000]
            True
            
            sage: mu = {0:[0,1], 1:[0]}
            sage: f = words.FixedPointOfMorphism(mu,0); f
            Fixed point beginning with 0 of the morphism {0: [0, 1], 1: [0]}
            sage: F = words.FibonacciWord(); F
            Fibonacci word over [0, 1], defined recursively
            sage: f[:10000] == F[:10000]
            True
        """
        # Some type checking.
        if first_letter != morphism[first_letter][0]:
            raise TypeError, "first_letter(=%r) must be the first element of morphism[first_letter](=%r)"%(first_letter,morphism[first_letter])
        if morphism[first_letter] == [first_letter]:
            raise TypeError, "first_letter(=%r) must be a proper prefix of morphism[first_letter](=%r)"%(first_letter,morphism[first_letter])
        # Define an iterator for the letters.
        def LetterIterator(morphism, first_letter):
            r"""
            Returns an iterator of the letters of the infinite word that starts
            with first_letter and is the fixed point of the morphism.

            If w is the word, then this iterator:
                outputs the elements of morphism[ w[i] ],
                appends morphism[ w[i+1] ] to w,
                increments i.
            """
            w = morphism[first_letter][:]
            i = 0
            while w:
                for a in morphism[w[i]]:
                    yield a
                else:
                    w += morphism[w[i+1]]
                    i += 1
            ## This appears to be slower than the above, and I'm not sure why.
            #while w:
            #    for a in morphism[w.pop(0)]:
            #        yield a
            #    w += morphism[w[0]]
        # Construct and return the word.
        w = w_0_1(LetterIterator(morphism, first_letter), format='iterator')
        w._morphism = morphism
        w.rename("Fixed point beginning with %r of the morphism %s"%(first_letter,morphism))
        return w

    def CodingOfRotationWord(self, alpha, beta, x=0):
        r"""
        Return the infinite word obtained from the coding of rotation of
        parameter $(\alpha,\beta, x)$.
        
        "The coding of rotation corresponding to the parameters $(\alpha,
        \beta, x)$ is the symbolic sequence $u = (u_n)_{n\geq 0}$ defined over
        the binary alphabet $\{0, 1\}$ by $u_n = 1$ if $x+n\alpha\in[0, \beta[$ 
        and $u_n = 0$ otherwise." (see[1])
        
        EXAMPLES:
            sage: alpha = 0.45
            sage: beta = 0.48
            sage: w = words.CodingOfRotationWord(0.45, 0.48); w
            Coding of rotation with parameters (0.450000000000000, 0.480000000000000, 0)
            sage: w[:100]
            Finite word of length 100 over [0, 1]
            sage: w[:30]
            word: 110101010100101010101101010101
        
        REFERENCES:
            [1] B. Adamczewski and J. Cassaigne, On the transcendence of real
            numbers with a regular expansion. J. Number Theory 103 (2003), pp.
            27-37.
        """
        def f(n):
            hauteur = x + n * alpha
            fracH = hauteur.frac()
            if fracH < 0:
                fracH += 1
            if 0 <= fracH < beta:
                return 1
            else:
                return 0
        w = w_0_1(f, format='function')
        w.rename("Coding of rotation with parameters (%s, %s, %s)"%(alpha,beta,x))
        return w
    
    def CharacteristicSturmianWord(self, cf, repeat=True):
        r"""
        Returns the characteristic Sturmian word with slope given by the
        continued fraction [0, cf[1], cf[2], ...]. Here cf should be an
        iterator.

        Definition. The characteristic Sturmian word of slope 
            alpha = [0, d[1] + 1, d[2], d[3], ...] 
        is the limit of the sequence 
            $s_0 = 1$, $s_1 = 0$, ..., $s_{n+1} = s_n^{d_n} s_{n-1}$ for $n > 0$.

        Equivalently, the n-th term of the characteristic Sturmian word of
        slope a is floor(a (n+1)) - floor(an) + floor(a).

        EXAMPLES:
            sage: def cf():
            ...     yield 0
            ...     while True: yield 1
            ...
            sage: F = words.CharacteristicSturmianWord(cf()); F
            Characteristic Sturmian word over [0, 1], defined recursively
            sage: Fib = words.FibonacciWord(); Fib
            Fibonacci word over [0, 1], defined recursively
            sage: F[:10000] == Fib[:10000]
            True
        """
        if not repeat:
            d = iter(cf)
        else:
            d = cycle(cf)
        def LetterIterator():
            s0 = [0]
            s1 = [1]
            s1, s0 = s1*(d.next()-1) + s0, s1
            n = 0
            while True:
                for i in s1[n:]:
                    n += 1
                    yield i
                else:
                    s1, s0 = s1*d.next() + s0, s1
        w = w_0_1(LetterIterator(), format='iterator')
        w.rename("Characteristic Sturmian word over %s, defined recursively" % w.alphabet().string_rep())
        return w
    
    def StandardEpisturmianWord(self, directive_word):
        r"""
        Returns the standard episturmian word (or epistandard word) directed by 
        the given "directive word". Over a 2-letter alphabet, this function gives 
        characteristic Sturmian words.
               
        INPUT: 
            directive word -- an infinite word or a period of a periodic infinite word         
        
        EXAMPLES:
            sage: Fibonacci = words.StandardEpisturmianWord(Word('ab')); Fibonacci
            Standard episturmian word over ['a', 'b']
            sage: Fibonacci[:25]
            word: abaababaabaababaababaabaa
            sage: Tribonacci = words.StandardEpisturmianWord(Word('abc')); Tribonacci
<<<<<<< .working
            Standard episturmian word over ['a', 'c', 'b']
=======
            Standard episturmian word over ['a', 'b', 'c']
>>>>>>> .merge-right.r169
            sage: Tribonacci[:25]
            word: abacabaabacababacabaabaca
            sage: S = words.StandardEpisturmianWord(Word('aabcabada')); S
<<<<<<< .working
            Standard episturmian word over ['a', 'c', 'b', 'd']
=======
            Standard episturmian word over ['a', 'b', 'c', 'd']
>>>>>>> .merge-right.r169
            sage: print S[:75]
            word: aabaacaabaaabaacaabaabaacaabaaabaacaabaaabaacaabaabaacaabaaabaacaabaadaabaa
            sage: S = words.StandardEpisturmianWord(Fibonacci); S
            Standard episturmian word over ['a', 'b']
            sage: S[:25]
            word: abaabaababaabaabaababaaba
            sage: S = words.StandardEpisturmianWord(Tribonacci); S
<<<<<<< .working
            Standard episturmian word over ['a', 'c', 'b']
=======
            Standard episturmian word over ['a', 'b', 'c']
>>>>>>> .merge-right.r169
            sage: S[:25]
            word: abaabacabaabaabacabaababa
            sage: words.StandardEpisturmianWord(123)
            Traceback (most recent call last):
            ...
            TypeError: directive_word is not a word, so it cannot be used to build an episturmian word
            sage: words.StandardEpisturmianWord(Words('ab'))
            Traceback (most recent call last):
            ...
            TypeError: directive_word is not a word, so it cannot be used to build an episturmian word

        """
        if not is_Word(directive_word):
           raise TypeError, "directive_word is not a word, so it cannot be used to build an episturmian word"
        elif is_FiniteWord(directive_word):
           d = cycle(directive_word)
        else:
           d = iter(directive_word)
           
        W = directive_word.parent()
               
        def LetterIterator():
            w = W(d.next())
            n = 0
            while True:
                  for x in w[n:]:
                      n += 1
                      yield x
                  else:
                      w = W(w*W(d.next())).palindromic_closure()
        
        epistandard = W(LetterIterator(), format='iterator')
        epistandard.rename("Standard episturmian word over %s" % epistandard.alphabet().string_rep())
        return epistandard
        
    def MinimalSmoothPrefix(self, n):
        r"""
        This function finds and returns the minimal smooth prefix of length n.
        
        INPUT:
            n -- the desired length of the prefix
            
        OUTPUT:
            word -- the prefix
            
        Be patient, this function can take a really long time if asked 
        for a large prefix.
        
        EXAMPLES:
            sage: words.MinimalSmoothPrefix(10)
            word: 1212212112
        """
        tab = []
        W = Words([1, 2])
        suff1 = W([1, 2, 2]).phi_inv()
        suff2 = W([2, 2]).phi_inv()
        w = [1]
        tab = _build_tab(1, tab, W)
        for k in xrange(1, n):
            if suff1._phi_inv_tab(tab) < suff2._phi_inv_tab(tab):
                w.append(1)
                tab = _build_tab(1, tab, W)
            else:
                w.append(2)
                tab = _build_tab(2, tab, W)
        return W(w)
    
    def RandomWord(self, n, m=2):
        """
        Returns a random word of length n on the alphabet {0,1,...,m-1}.
        
        INPUT:
            n -- integer, the length of the word
            m -- integer (default 2), the size of the output alphabet
        
        EXAMPLES:
            sage: words.RandomWord(10)         # random results
            word: 0110100101
            sage: words.RandomWord(10, 4)      # random results
            word: 0322313320
        """
        return Words(xrange(m))([randint(0,m-1) for i in xrange(n)])

words = WordGenerator()

class ChristoffelWord_Lower(FiniteWord):
    r"""
    Returns the lower Christoffel word of slope p/q, where p and q are
    relatively prime nonnegative integers.
    
    
    The Christoffel word of slope p/q is obtained from the Cayley graph of
    Z/(p+q)Z with generator q as follows. If u --> v is an edge in the
    Cayley graph, then v = u+p mod (p+q). Label the edge u --> v by 1 if u
    < v and 0 otherwise. The Christoffel word is the word obtained by
    reading the edge labels along the cycle beginning from 0.
    
    EXAMPLES:
        sage: from sage.words.word_generators import ChristoffelWord_Lower
        sage: w = ChristoffelWord_Lower(1,0); w
        Lower Christoffel word of slope 1/0 over the alphabet [0, 1]
        sage: print w
        word: 1
        
        sage: w = ChristoffelWord_Lower(0,1); w
        Lower Christoffel word of slope 0/1 over the alphabet [0, 1]
        sage: print w
        word: 0
        
        sage: w = ChristoffelWord_Lower(1,1); w
        Lower Christoffel word of slope 1/1 over the alphabet [0, 1]
        sage: print w
        word: 01
        
        sage: w = ChristoffelWord_Lower(4,7); w
        Lower Christoffel word of slope 4/7 over the alphabet [0, 1]
        sage: print w
        word: 00100100101
    """
    def __init__(self, p, q):
        # Compute gcd of p, q; raise TypeError if not 1.        a, b = p, q        while a:            a, b = b%a, a        if b != 1:            raise TypeError, "%s and %s are not relatively prime"%(p, q)
        # Compute the Christoffel word
        w = []
        u = 0
        if (p, q) == (0, 1):
            w = [0]
        else:
            for i in range(p + q):
                v = (u+p) % (p+q)
                new_letter = 0 if u < v else 1
                w.append(new_letter)
                u = v
        super(ChristoffelWord_Lower, self).__init__(
                            MakeContent(w, id_f, format='list'), w_0_1)
        self.__p = p
        self.__q = q
        self.rename("Lower Christoffel word of slope %s/%s over the alphabet %s" % (self.__p, self.__q, self.alphabet().string_rep()))

    def markoff_number(self):
        r"""
        Returns the Markoff number associated to the Christoffel word self.
        
        The Markoff number of a Christoffel word w is 1/3 trace M(w), 
        where M(w) is the 2x2 matrix obtained by applying the morphism
        
            0 -> matrix(2,[2,1,1,1])
            1 -> matrix(2,[5,2,2,1])
        
        EXAMPLES:
            sage: from sage.words.word_generators import ChristoffelWord_Lower
            sage: w0 = ChristoffelWord_Lower(4,7) 
            sage: w1, w2 = w0.standard_factorization()
            sage: (m0,m1,m2) = (w.markoff_number() for w in (w0,w1,w2))
            sage: (m0,m1,m2)
            (294685, 13, 7561)
            sage: m0**2 + m1**2 + m2**2 == 3*m0*m1*m2
            True
        """
        from sage.matrix.constructor import matrix
        eta = {0:matrix(2,[2,1,1,1]), 1:matrix(2,[5,2,2,1])}
        return self.apply_morphism(eta).trace()/3
    
    def standard_factorization(self):
        r"""
        Returns the standard factorization of self.
        
        The standard factorization of a Christoffel word w is the unique
        factorization of w into two Christoffel words.
        
        EXAMPLES:
            sage: from sage.words.word_generators import ChristoffelWord_Lower
            sage: w = ChristoffelWord_Lower(5,9) 
            sage: print w
            word: 00100100100101
            sage: w1, w2 = w.standard_factorization()
            sage: print w1
            word: 001
            sage: print w2
            word: 00100100101
            
            sage: w = ChristoffelWord_Lower(51,37) 
            sage: w1, w2 = w.standard_factorization()
            sage: w1
            Lower Christoffel word of slope 11/8 over the alphabet [0, 1]
            sage: w2
            Lower Christoffel word of slope 40/29 over the alphabet [0, 1]
            sage: w1 * w2 == w
            True
        """
        p, q = self.__p, self.__q
        index = 0
        u = 0
        for i in range(p + q):
            v = (u+p) % (p+q)
            if v == 1:
                index = i
                break
            u = v
        w1, w2 = self[:index+1], self[index+1:]
        return Factorization([ChristoffelWord_Lower(w1.count(1),w1.count(0)),
                ChristoffelWord_Lower(w2.count(1),w2.count(0))])
